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1.  INTRODUCTION 

Elliptic  and  hyperbolic  partial  differential  equations  are,  at  the  heart  of  most 
mathematical  models  used  in  engineering  and  physics,  giving  rise  to  extensive 
computations.  Often  the  problems  that  one  would  like  to  solve  exceed  the 
capacity  of  even  the  most  powerful  computers.  On  the  other  hand,  the  time 
required  is  too  great  to  allow  inclusion  of  advanced  mathematical  models  in  the 
design  process.  Also  iterative  processes  for  solving  the  algebraic  equations 
arising  from  discretizing  partial-differential  equations  are  stalling  numerical 
processes,  in  which  the  error  has  relatively  small  changes  from  one  iteration  to 
the  next.  The  computer  grinds  very  hard  for  very  small  or  slow  real  physical 
effect  with  the  use  of  too-fine  discretization  grids.  In  this  case,  in  large  parts  of 
the  computational  domain,  the  mesh  size  is  much  smaller  than  the  real  scale  of 
solution  changes.  In  particular,  convergence  of  iterative  methods  for  elliptic  grid 
generation  based  on  non-linear  grid  generation  equations  is  extremely  slow. 

For  the  above  reason,  the  authors  have  worked  on  the  improvement  of  the 
non-linear  elliptic  grid  generation.  In  this  paper,  the  Laplace  equation  and  the 
algebraic  transformation  were  presented  for  the  domains  in  2D  and  3D  physical 
space.  The  second  order  finite  difference  scheme  was  used  for  the  discretization 
of  the  grid  generating  equations.  The  linear  system  is  solved  by  the  ADI  method 
and  the  convergence  was  accelerated  by  the  multigrid  FAS  scheme.  The 
performance  characteristic  of  the  algorithm  was  discussed  and  illustrations  were 
made. 

2.  2D  AND  3D  ELLIPTIC  GRID  GENERATION 

The  elliptic  grid  generation  method  used  in  this  paper  is  based  on  the  use  of  a 
composite  mapping  [1],  It  is  a composition  of  an  algebraic  transformation  and  an 
elliptical  transformation  based  on  Laplace  equations.  The  algebraic 
transformation  is  a differential  one-to-one  mapping  from  computational  space 
onto  a parameter  space.  The  parameter  space  and  the  computational  space  are 
unit  squares.  The  algebraic  transformation  will  only  depend  on  the  prescribed 
boundary  grid  point  distribution.  The  control  functions  are  defined  based  on  the 
algebraic  transformation.  The  elliptical  transformation  is  a differential  one-to-one 
mapping  from  parameter  space  onto  the  physical  domain.  The  elliptical 
transformation  depends  only  on  the  shape  of  the  domain  and  is  independent  of 
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the  prescribed  boundary  grid  point  distribution.  The  composition  of  these  two 
mappings  defines  the  interior  grid  point  distribution  and  is  a differentiable  and 
one-to-one  for  2D  domains  and  surfaces  and,  probably,  also  for  3D  domains. 

2.1  Grid  Generation  Equations 

For  2D  problem,  the  grid  generating  system  of  elliptic  partial  differential 
equations  are  as  follow: 

Pxg  + 2 QXfr  + Rxm  + SXf  + Txn  = 0 (1 ) 

Where 

P = (xn , xn ),Q  = -(*£  ,xv),R  = ( x ( , ), 

S = PPt\  + 2<2P,2  + RP22 , T = PP*  + 2 QPtl  + RP>2  (2) 

The  control  functions  are  given  by 
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and  the  matrix  T is  defined  as 
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The  six  coefficients  of  the  vectors 

P\  1 = (^ii  > P\i  Y » P\2  — (^12  > Pn  Y andP22  — (P22  ,P22)t  are  called  the  control 
functions.  The  two  algebraic  equations  that  define  the  transformation  are  given 

by 

s = sEi(Z)(l-t)  + sEi(%)t  (4) 


and  t = tEx  (rj)(  1 -s)  + tE2  ( Tj)s 


(5) 


The  3D  grid  generating  system  is  too  complicated  to  describe  here.  It  can  be 
found  in  [1]. 

2.2  Discretization  Method 

Let’s  use  2D  case  as  an  example  and  consider  a uniform  rectangular  grid  of 
size  (A  + 1)(M+1)  defined  as  ^ij=^i=i/N,  Tji  . =T]j  = j / M, 
i = 0...N\j  = 0...M. 
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Let  X i . be  prescribed  on  the  boundary  of  this  grid.  We  are  going  to  compute 
X i . in  the  interior  of  the  computational  grid  based  on  the  solution  of  the 
Poisson  system  defined  by  Eq.  1 . 

The  solution  of  this  system  of  nonlinear  elliptical  equations  is  obtained  by 
Picard  iteration. 


/>k-1xk^  + 20k-‘xk^  + Rk-1xk 
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(xk\  Xk  ln),  Qkl  = -(xk\  xk\),  Rk  l = (xk\  xk~\). 


Sk'1  = Pk']  p\i  + 2QkA  p\2  + RkA  P'22, 14"'  = Pl 


'k-'  P2n  + 2Qk-'  P2l2  + Rk-'  P222. 


The  six  control  functions  P\\,  P\2,  P'22,  P2n,  P2 \2,  P2 22  are  computed 
according  to  the  equations  given  by  Eq.2  and  by  applying  second  order  central 
difference  schemes  for  the  discretizations  of  s^,s  ,s*  ,t„,tnn, 

> s£,  'Sn'h'*ri 


The  arc  length  normalized  variables  (st  j , tt  j ) at  the  boundary  are  computed 
as  follows: 

• Compute  the  distance  (/)  between  succeeding  points  at  the  each  boundary. 

• Define  the  length  of  the  edge  (L)  along  each  boundary  as  the  sum  of  the 
distances  between  succeeding  points  along  that  boundary. 

• The  normalized  distance  ( d ) between  succeeding  points  is  then  given  by  HL. 

• The  arc  length  normalized  variables  s,  ; and  tt  j at  the  boundary  are  defined 
by 

s0j=  0,  sNJ=l,  j = 0...M , ti0=  0,  tiM  = 1,  i = 0...N, 
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The  arc  length  normalized  variables  (5;  • , tt  j ) in  the  interior  of  the  grid  are 
computed  according  to  the  algebraic  straight-line  transformation  given  by  Eq.(3) 
and  (4).  Simultaneously  solving  the  two  linear  algebraic  equations  yields 


Sij  = J/,0  (1  - t,,j ) + s,,M  t,.j . h.j  = *o.j  (!  - si,j ) + 1. 


N.jSiJ ’ 


at  each  node  (i,  j ) e (1  ...N  - 1;1...M  - 1) 

The  steps  for  an  iteration  to  improve  the  current  approximation  xk~'  is  as 
follows: 

• The  coefficients  Pk~'  ,Qk  ' , Rk  ' ,Sk~l  ,Tk~'  are  computed  by  applying  central 
difference  schemes  for  the  discretization  of  x'T'1  and  xk~'  . The  six  control 
functions  remain  unchanged  during  the  iterative  process. 
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• Using  central  difference  schemes  to  discretize  , x*v , xji , x* . The 

discretization  of  the  mixed  derivative  using  central  difference  schemes  is  done  in 
a way  described  in  [2]. 

• We  then  obtain  a linear  system  of  equations  for  the  unknowns 
X-j,i  = 0 ...N;  j = 0...M  with  Dirichlet  boundary  conditions.  A Multigrid  solver 
is  used  to  solve  this  linear  system.  The  solver  is  called  twice  to  compute  the  two 
components  jc*  • and  y* . . 

The  initial  approximation  a:0  is  obtained  by  an  algebraic  grid  generation.  The 
above  iterative  process  is  repeated  until  a required  approximation  to  the  solution 
has  been  obtained. 

The  discretization  method  for  3D  problem  can  also  be  found  in  [1]. 

3.  FINITE  DIFFERENCE  AND  FINITE  VOLUME  METHOD 


In  this  paper,  the  second  order  finite  difference  scheme  was  used  to  discretize 
the  equations  and  the  boundary  conditions.  For  2D  problem,  Eq.  (1)  becomes: 
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Grids  obtained  by  the  nonlinear  elliptic  Poisson  grid  generation  system 
defined  by  Eq.(l)  are  grid  folding  free  and  have  an  excellent  interior  gird  point 
spacing  distribution.  However,  the  computed  grids  are  in  general  not  orthogonal 
at  the  boundary,  but  we  need  grids  to  be  orthogonal  at  the  boundary,  especially 
for  Navier-Stokes  computations.  The  orthogonality  of  the  gird  in  a boundary 
layer  is  often  desired. 

The  s coordinate  in  parameter  space  P satisfies  the  linear  second-order 
elliptic  equation  (JallSg  + + (Jal2s^  + Ja22s  )n  = 0 . The  t 

coordinate  is  obtained  in  the  same  way.  A finite-volume  cell-centered  nine-point 
stencil  approach  is  used  to  obtain  the  discretized  equations  for  boundary 
orthogonality  of  2D  problem  and  seven-point  stencil  for  that  of  3D  problem. 

4.  MULTIGRID  TECHNIQUE 


As  we  know,  Gauss-Seidel  relaxation  for  solving  Eq.6  typically  stalls  after  a 
few  iterations.  This  is  because  Gauss-Seidel,  though  effective  for  high-frequency 
errors  components,  has  very  little  effect  on  low-frequency  components.  Multigrid 
capitalizes  on  this  “smoothing”  property  of  Gauss-Seidel  by  visiting  coarser  grids 
to  resolve  smooth  errors.  But  the  regular  multigrid  correction  scheme  is  only 
valid  for  linear  problem.  In  order  to  accommodate  the  nonlinearities  in  Eq.6,  we 
applied  the  Full  Approximation  Scheme  (FAS)  version  of  multigrid  with  bilinear 
interpolation  and  full  weighting  of  residuals  and  approximations. 
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4.1.  Full  Approximation  Scheme  (FAS) 

In  non-linear  problem,  the  difference  Lhuh  - LhUhcm  no  longer  be  replaced 
by  Lhvh , where  vh  is  the  truncation  error  onfi|)e_grid.  Hence,  we  ii^pduce  a new 
coarse  grid  variable  U21,  defined  as  u n,  = 1 1,  ui,  +v2/l,  where  1 1,  represents  a 
so-called  restriction  operator  which  interpolates  fine  grid  solution  variables  to  the 
coarse  grid  and  v2/l  is  the  truncation  error2;jon  coarse  grid.  The  coarse  grid 
equation  can  be  written  as  L2hui>,  - Llhl h Ui,  - —lh ' rh  , lh ' represents  the 
restriction  operator  which  transfers  residua^  from  fine  to  coarse  grids  and  rh  is 
the  residue  on  fine  grid.  The  operators  lh  may  in  principle  be  different  from 
each  other. 

It  is  useful  to  rewrite  the  above  equation  as: 

— — 2 /i  — 2/j 

.Z^2/j ^ 2/1  — S 2)1 » where  S 2^  l^2h ^ ^ Hu  I ^ . 

In  this  form,  the  coarse  gird  equation  is  seen  to  take  on  a similar  structure  to  the 
original  fine  grid  equation,  with  a modified  source  term.  This  enables  us  to  use 
similar  techniques  for  solving  both  the  coarse  and  fine  grid  problems.  Once  the 
coarse  grid  equations  have  been  solved,  either  exactly  or  approximately,  the  fine 
grid  variables  are  updated  as: 

— new  — old  . —new  —2h—old 

Uh  — li  h +/2/,(m  2/i  —lh  Uh  ) 

— new  — old  ^ ^ 

which  can  also  be  written  as  Ui,  =Uh  + I2hv2h  ■ 1 2 1,  *s  calUd  interpolation 
operator. 

The  FAS  method  is  carried  out  in  the  following  steps.  [4] 

• Restrict  the  current  approximation  and  its  fine-grid  residual  to  the  coarse  grid: 
r2h  =Ihh(fh  ~ Lh  (vi, ))  and  V2i,  = ll’'vh  > where  fh is  the  right-hand  term  of 
nonlinear  equation  system  on  fine  grid. 

• Solve  the  coarse-grid  problem  Llh  ( u2h ) = L2h  (v2h  ) + r2h. 

• Compute  the  coarse-grid  approximation  to  the  error:  e2h  =u2h  ~v2h. 

• Interpolate  the  error  approximation  up  to  the  fine  grid  and  correct  the  current 
fine-grid  approximation:  vh  vh  + I2he2h , where  h:  the  mash  space  on  finest 
grid,  2h:  the  mesh  space  on  coarse  grid. 

4.2.  The  Solver  for  Linear  Systems 

The  linear  systems  were  solved  by  applying  alternating  direction  implicit 
methods  (ADI).  The  basic  idea  for  ADI  method  is  as  follows:  a complete 
iteration  consists  of  two  scans.  The  first  scan  goes  from  the  smallest  index 
number  to  the  largest  one  from  one  row  to  the  other  followed  by  the  same  in  the 
column  direction.  The  second  scan  is  almost  the  same  as  the  first  one,  except  for 
going  from  the  largest  index  number  to  the  smallest  one.  This  method  can 
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increase  the  convergence  rate  at  the  cost  of  some  complication  in  the 
computational  algorithm.  [3] 

To  illustrate  the  advantage  of  FAS  method,  the  author  compared  ADI  method 
with  FAS  and  ADI  method. 

5.  ILLUSTRATIONS 
5.1.  For  2-D  Case 

In  Fig.  2(a),  the  body-fitted  C-grid  around  a NACA  0012  airfoil  is  displayed. 
The  close-up  of  mesh  near  the  airfoil  is  shown  in  Fig.  2(b)  where  grids  are 
orthogonal  at  the  boundaries.  The  mesh  is  clustered  near  the  wall  in  the  wall- 
normal  direction  ( Tj ) as  well  as  in  the  vicinity  of  the  leading  and  trailing  edge  of 
the  airfoil  in  the  direction  parallel  to  the  wall  (£) . The  grid  contains  769  points 
in  the  £ direction  and  129  points  in  the  T]  direction. 

As  seen  in  the  figures,  the  grid  generated  is  grid  folding  free  and  the  interior 
grid  point  distribution  is  a good  reflection  of  the  prescribed  boundary  grid  point 
distribution.  The  initial  grid  is  obtained  with  algebraic  grid  generation  and  is 
required  as  the  initial  solution  for  the  non-linear  elliptical  Poisson  system.  The 
final  grid  is  independent  of  the  initial  grid.  The  quality  of  the  initial  grid  is 
unimportant  and  severe  grid  folding  of  the  initial  grid  is  allowed. 

Fig.  3.  shows  the  Log  (residue)  vs.  Work  Units  of  the  ADI  and  the  multigrid 
method  for  2D  case.  From  the  figure,  the  multigrid  method  reduced  the  residue 
significantly,  which  is  faster  than  the  ADI  method  for  the  above  problem.  The 
multigrid  method  saves  at  least  23  times  in  CPU  times  (or  work  unit)  than  the 
ADI  method. 

5.1.  For  3-D  Case 

The  geometry  of  the  delta  wind,  taken  from  the  experimental  work  of  Rieley 
& Lowson  (1998),  is  shown  in  Fig.  4.  The  sweep  angle  denoted  by  A is  85°  and 
the  leading-edge  angle  denoted  by  cr  is  30° . The  chord  length  is  taken  as  the 
characteristic  length  L , such  that  the  non-dimensional  chord  length  is  c = 1 ,0L . 
The  non-dimensional  thickness  of  the  delta  wing  is  h — 0.024L . 

In  Fig.5,  the  grids  around  the  delta  wind  are  being  shown.  The  mesh  is  H-type 
in  the  meridian  section  and  C-type  in  the  cross  section.  The  grids  are  orthogonal 
on  the  delta  wing  surface.  The  meshes  are  1 29  x 65  x 65 , where  the  sequence  of 
numbers  is  corresponding  to  the  axial,  the  spanwise  and  the  wall-normal 
direction,  respectively. 

Fig.  6.  shows  the  Log  (residue)  vs.  Work  Units  of  the  ADI  and  the  multigrid 
method  for  3D  case.  Great  improvement  of  convergence  speed  can  also  been 
seen  from  this  figure,  although  it  is  not  as  apparent  as  that  in  2D  case.  The 
multigrid  method  saves  at  least  2 times  in  CPU  times  (or  work  unit)  than  ADI 
method. 

6.  CONCLUSIONS 
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The  Multigrid  method  is  applied  to  elliptic  partial  differential  equations.  As 
seen  from  the  illustrations,  convergence  of  the  2D  and  3D  elliptic  grid  generation 
problems  is  greatly  improved  by  using  the  multigrid  solver,  although 
convergence  of  three  dimensional  elliptic  grid  generation  is  still  not  fast  enough 
and  need  to  be  improved.  Second  order  central  difference  schemes  used  for 
discretizing  the  grid  generation  equations  are  close  to  the  exact  differentiation 
and  provide  better  grid  point  distribution. 
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Computational  apace  C Parameter  space  P Domain  D 

FIG.  1.  Transformation  from  computational  (£,77)  space  to  a domain  D in 
Cartesian  (A,  y)  space. 
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Comparison  of  Log  (residue)  vs.  Work  Units 
of  ADI  and  multigrid  method 


Work  Units 


FIG.  3.  Comparison  of  Log  (residue)  vs.  FIG.  4.  Schematic  of  the  delta  wing 
Work  Units  of  ADI  and  multigrid 
method  in  2D  case 


FIG.  5.  H-C  type  grid  around  a 85°  FIG.  6.  Comparison  of  Log(residue)  vs. 

sweep  delta  wing  Work  Units  of  ADI  and 
multigrid  method  in  3D  case 


